Statics and dynamics of BEC's in double square well potentials. 
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^5 , In this paper we treat the behavior of Bose Einstein condensates in double square well potentials, 

f ^ ■ both of equal and different depths. For even depth, symmetry preserving solutions to the relevant 

' nonlinear Schrodinger equation is known, just as in the linear limit. When the nonlinearity is strong 

enough, symmetry breaking solutions also exist, side by side with the symmetric one. Interestingly, 
solutions almost entirely localized in one of the wells are known as an extreme case. Here we outline 
a method for obtaining all these solutions for repulsive interactions. The bifurcation point at which, 
for critical nonlinearity, the asymmetric solutions branch off from the symmetry preserving ones 
is found analytically. We also find this bifurcation point and treat the solutions generally via a 

in . 

Josephson Junction model. 

When the confining potential is in the form of two wells of different depth, interesting new 
phenomena appear. This is true of both the occurrence of the bifurcation point for the static 
' solutions, and also of the dynamics of phase and amplitude varying solutions. Again a generalization 

of the Josephson model proves useful. The stability of solutions is treated briefly. 
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I. INTRODUCTION 



The nonlinear Schrodinger equation is a powerful tool for describing Bose Einstein condensates at zero temperature. 
Double well potentials are an important class of configurations to which this tool can be applied. For square wells, 
exact solutions are known to exist 0. In Zin et al. 0, we outlined a method for obtaining such exact solutions for 
a symmetric double square well situation with attractive interaction. We found an exact criterion to determine the 
bifurcation point. Here we perform similar calculations for the repulsive case and extend our treatment, including 
wells of different depth but also a stability analysis. The repulsive case is perhaps more interesting in view of the fact 
I ' that, for this interaction, situations such that most of the condensate was contained in one of the wells have been 
seen experimentally 0- We also present some dynamic calculations not included in Q for both kinds of interaction. 
t^J- ' These lean somewhat on a generalized Josephson Junction model 0, IE IE 13 • We note in passing that a Josephson 
r*** ! Junction for a Bose-Einstein condensate was first obtained by Inguscio's group 0, see also 

Symmetry breaking solutions that are known to exist for positive nonlinearity (repulsive interaction in the case 
of BEC, dark solitons in nonlinear optical media) often tend to localize the wave function in one of the wells. This 
happens for the nonlinearity exceeding a critical value at which the asymmetric solutions branch off from the symmetry 
preserving ones in the parameter space. Therefore, we can talk about bifurcation at this critical value of nonlinearity. 
The existence of such solutions of the nonlinear Schrodinger equation was first pointed out in the context of molecular 
states for repulsive interaction |lOj| . as will be treated here. Importantly, the effect of this spontaneous symmetry 
breaking has been observed in photonic lattices . It should be stressed that the nature of bifurcation depends on 
I 1 the symmetry of the problem and is of the pitchfork variety for even wells and saddle point for uneven wells |T^ . 

In this paper we will consider a double square well potential, first symmetric and then asymmetric. The asymmetric 
potential leads to more complicated profiles. As far as we know, these square well configurations are the only ones 
for which exact solutions exist. These solutions are all in the form of Jacobi elliptic functions. One of the problems 
considered here, extending to repulsive interaction, is how to proceed from easily obtainable symmetric double well 
solutions of the linear Schrodinger equation to the fully nonlinear case, and from so obtained symmetric solutions on 
to the bifurcated, asymmetric ones. When the potential is asymmetric both bifurcation of the static solutions and the 
dynamics of oscillating solutions will be seen to become very different from those for the symmetric potential. 

The manuscript is composed as follows: In section^Jwe derive symmetry preserving states starting from the linear 
limit and then gradually increasing the nonlinear interaction. In section III we investigate the symmetry breaking 
states that branch off from the symmetry preserving ones in the parameter space. We give a simple exact formula 
for the bifurcation point. Section IV treats asymmetric potentials. Section V is devoted to dynamics treated by the 
Josephson model, particularly useful at the bifurcation point, and then numerically. Results are consistent by all 
three methods (sections III, IV and V). Some concluding remarks wind up the text (section VI). Heavier calculations 
have been relegated to the Appendix. 

This paper can be read independently of reference . 
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II. ANTISYMMETRIC STATES FROM THE LINEAR LIMIT (SYMMETRIC WELLS) 



We start from the nonlinear Schrodinger equation 



s 
dx 2 



^-+V(x)+r,\f(x)\ 2 



/0) = nf(x) 



Here the potential is of the form 



V for \x\ < b 
V(x) = { for b < \x\ < a 
oo for \x\ > a 



(1) 



(2) 



See Fig. ^ Solutions in the three regions will be written as 



fi(x) for — a < x < —b 
f(x) = { f 2 (x) for \x\ < b 

fs{x) for b < x < a 



(3) 



The solutions vanish on and outside the outer boundaries |a;| > a. We assume continuity of f[x) and its derivative at 
x = ±6 and normalization to J^ a \f(x)\ 2 dx = 1. The symmetric solutions are 



fi(x) — A sn(k(x + a)\m) 
f 2 {x) = A 2 nc(k 2 x\m 2 ) 
fz(x) — — A sn(k(x — a)\m) 

and the antisymmetric solutions, which will be of particular interest here, are 

fi(x) — A sn(k(x + a)\m) 
h{®) = ~ A 2 sc{k 2 x\m 2 ) 
fs[x) — A six{k(x — a)\m) . 



(4) 



(5) 



Here f\{x) and fs{x) have been chosen to be zero at the ends, and also so as to preserve even and odd parity 
respectively for the two cases. The parameters of the symmetric solutions are found from Eq. to satisfy (Vq > fi) 



A 1 



2mk 2 



A2 _ 2(1 - m 2 )k 2 2 



and for the antisymmetric solutions we have 



A = 



2mk 2 
V 



2(1 - m 2 )k\ 



// = (! + m)k 2 = (1 - 2m 2 )k 2 + V 



H = (1 + m)k 2 = (m 2 - 2)kl + V . 



(6) 



(7) 



Positive roots for all the A's are taken throughout. We choose /J,, m and m 2 to generate all the other constants. These 
three parameters will determine the solution completely. 

We now concentrate on the antisymmetric case, as we have checked that bifurcation only occurs for this case in the 
lowest mode as suggested by Fig. 0](and also by Fig. 3 of Ref. [13j)- We have two continuity conditions at ±6 



Asn(koj\m) = —A 2 sc(—k 2 b\m 2 ) 

Akcn{kuj\rn)dn(kuj\m) = — A 2 k 2 dc(— k 2 b\m 2 )nc(— k 2 b\m 2 ). 
Here u> — a — b. The normalization of the wave function, J" a da; |/(x)| 2 = 1 yields 



(8) 
(9) 



2A 2 / sn 2 (k(x — a )\m)dx + 2Al / sc 2 (fc 2 

J a JQ 



The above normalization condition works out as: 

4fc 2 [sn(fc 2 6|m 2 ) dc(fc 2 6|m 2 ) - E(k 2 b\m 2 )} + Ak [ku - E(ku\m)} = 77, 



(10) 
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FIG. 1: Symmetric double square well potential. In the following 2a — 1, 2b — 0.1 and Vo = 300. 



where E(u\m) is the elliptic function of the second kind fLU ]. We now have three equations for m, m 2 and fi as 
required. Here a, b, Vq and n are fixed and describe one specific experimental setup. Up to now, the shooting method 
was used to find solutions pj. 

To systematically solve our equations we first turn to the linear limit r) — 0, k 2 = /x, k\ = Vo — \i. The functions 
are now easy to calculate 

fi(x) = A sinfc(a; + a) 
f2{x) = —A 2 sinhfc 2 a; 
fz{x) = A sinfc(a; — a). 



The two continuity conditions are: 



A sin kuj = A 2 sinhfc 2 & (H) 
Ak cos few = —A 2 k 2 coshfc 2 6, , (12) 



This last equation will give us linear /x in terms of fixed parameters. The normalization condition is 

sin2fco; sin ku f sinh2kib \ \ 
2k + sinh 2 k 2 b V + 2fc 



^ 2 -(--^ + ^T(-^+=^)) • (13) 



Giving a value for A, and A 2 follows from continuity. We are now ready to tie this solution up to the small r; limit in a 
perturbative manner. We notice that A, A 2l k, k 2 obtained in the linear approach become a zero order approximation 
in an r\ expansion. The parameters m and 1 — m 2 arc both of order r\ and follow from equations © and (|llfl 



A 2 sin 2 kuj A 2 

Linear /i, denoted by /xq, is found as the lowest root of 



2k 2 sinh 2 k 2 b2k 2 ' K J 



fi cot(,/jlb~u) + \/Vo - /x coth(vT^) - fi a b) = 0. (15) 



There will also be a small rj correction A^i such that /x = /io + A/x. This will complete the calculation of the three 
unknowns /x, m, m 2 in the small rj limit (Appendix). 

Now that we have a starting point, we can generate all symmetric solutions by gradually increasing r\. We introduce 
the notation /x = toq, m = mi, and m 2 . We write the conditions (jSJ) and 10 in the symbolic functional form 



h (m , mi,m 2 ) = k^Jm\ sn(fca;|mi) — k 2 ^/(l — m 2 ) sc(fc 2 &|m 2 ) = 0, (16) 
hi(mo, mi, m 2 ) — k 2 ^Jm\ cn(kui\mi) dn(/cw|mi) + k\\J\ — ra 2 dc(fc 2 6|m 2 ) nc(fc 2 6|m 2 ) = 0. (17) 

Here we used Eq. (JHJ to express the amplitudes A and A 2 in terms of the mi and the wavevectors k and k 2 , which 
in turn can be expressed in terms of the m^. The left hand side of equation (|10fl defines h 2 (mo,mi,m 2 ), which is 
evidently free of 77. Upon defining 770 = 0, 771 = 0, r\ 2 = r\ we write all three conditions ©, 10 and (|1L)|) simply as 

hi(mo,mi,m 2 ) = r/i for i = 0,1,2. (18) 
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In all three equations l |18| l functions ft,(mo, m±, m 2 ) on the left, remain free of 77. Hence if we increase r\ by a small 
increment A77 the parameters m, will increase by Am, governed by 

(££) - A "» (19) 

where Ar/i = (0,0, A17). Ass uming the matrix f g TO '. 1 to be nonsingular we can now generate increments in m^ by 
gradually increasing the control parameter 77. Inverting Eq. I|19|) we find 



III. SYMMETRY BREAKING STATES (SYMMETRIC WELLS) 

Even when the double well is symmetric, the nonlinear Schrodinger equation is known to admit symmetry breaking 
states. This is in contradistinction to the linear version, admitting only symmetric and antisymmetric states as treated 
in section ifTTjl. These symmetry breaking states are possible above a critical value of r\. They are of considerable 
physical interest, as they include situations such as the location of most of the wavefunction in one half of the double 
well. More generally, there is the possibility of very different profiles in the two halves. The solutions corresponding 
to symmetry breaking are known to bifurcate from the antisymmetric ones. Here we will give a condition defining the 
bifurcation points in parameter space and investigate how this bifurcation can be interpreted. We will give diagrams 
to illustrate this. Similar diagrams for a quartic potential can be found in however they do not correspond to 
any analytic solutions known to us. 

Solutions generalizing the symmetric case are: 

fi(x) = A\ sn{k\(x + a)\m{) 
f 2 (x) = A 2 nc(k 2 (x + d)\m 2 ) 
fa(x) = ~A 3 sn(k 3 (x - a)|m 3 ) 
and the generalization for the antisymmetric case is: 

fi(x) — A\ sn.{k\[x + a)\mi) 
h{ x ) = -A 2 sc(k 2 (x + d)\m 2 ) 
f 3 (x) = A 3 sn(k 3 (x - a)|m 3 ). 

When d = 0, toj = m 3 = m the solutions Q and J^J) are recovered. Once again we concentrate on a generalization 
of the antisymmetric case, as the only one branching off from a basic mode. 

We now have five conditions for [i, mi, m 2 , m 3 , d, which we will denote mj I — 0, ...4 and in place of equation Ijl8(l 
we have gj = rji, see the Appendix. One solution is d = O.mi = m 2 = m, as we know, and conditions JHJ, © an d 
(|10ll are recovered. However, as we will see, above a certain threshold in 77 a second solution appears. The value of 77 
at which this bifurcation occurs will be denoted by r\uf- The second solution branches off the antisymmetric one at 
this point. To find it we note that at such a point the antisymmetric solution is continuous with respect to 77, whereas 

the asymmetric one is not. Therefore we expect the 3x3 matrix (j^-^J to be nonsingular, whereas the 5x5 matrix 

(Smj) ^ e sm S mar & t this point. Simple algebra shows that the determinant of the 5x5 matrix can be factorized 
at the bifurcation point for which mi = 777,3 and d = 0. 

det ( pL) = det ( p.) D 2 (21) 



and D 2 is found to be given by: 



Do 



\dmj J \dm>3 

dgo dg 2 _ dg dg 2 
dmi dd dd dui\ 



(22) 

mi— m,3,d— 

In view of the above, D 2 = 0. This condition can be expressed in terms of variables characterizing the antisymmetric 
solution. If we write c 
the bifurcation point: 



solution. If we write conditions (jl6|l . I|17|l as ho = h { 1} - 4 2) and h i = h i 1] + h i ] 

we obtain a simple condition for 



dbh^dmh^ + d m h ( ^d b h { 2) = 0. (23) 
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FIG. 2: Antisymmetric and asymmetric solutions for the same value of 77 = 10. 



This further simplifies to: 



dm 



\mk A — TOfc 2 Vbsn 2 (fctj|m)] = 0. 



(24) 



Thus we can find the bifurcation point on the antisymmetric branch in terms of just two of the antisymmetric 
variables. Having that point we can move out onto the asymmetric branch using: 



Am/ = 



dgi 
dm , 



r\i = for 7 = 0,1,2,3 and 774 = 77. 



(25) 



This equation can be used everywhere except at the branch point, where second derivatives must come in. We have 
checked numerically that condition (|24|l is always satisfied at the bifurcation point. For an illustration of a bifurcated 
asymmetric solution paired with the corresponding symmetry preserving one far from bifurcation, see Fig. [21 

IV. ASYMMETRIC POTENTIAL 

Suppose now that the right hand well is somewhat shallower that on the left (V = V 3 , b < x < a, < V3 < /x). 
Otherwise we keep the notation of Fig. ^ Eqn. © is now replaced by 



2m\k\ 
V 



2(1 ™ 2)kl , A 2 = ^|M, n=(l + mi )kj = (l-2m 2 )k 2 2 + V = (l + m 3 )kl + V 3 . (26) 



Equation J7J) is similarly modified. When this is done all formulas for the symmetry breaking case formally carry 
through, with the understanding that fc| is now (/i — V3 ) / ( 1 + m 3 ) . Interchanging mj and m 3 no longer gives a trivial 
alteration. Solutions with most of the condensate on the left or on the right are no longer mirror images. Also, the 
linear limit is altered, relevant equations becoming 

fi(x) — A\ sinfci(a; + a) 
f2{x) = —A2 smhk 2 (x + d) 
f 3 (x) = A 3 sinfc 3 (a; - a). 

Even in this limit, we now have five equations for five unknowns: mi, m,2, m 3 , /x and d. This limit is thus nolonger 
much simpler than the fully nonlinear case. However this is not worth pursuing, as the interesting bifurcation does 
not now occur from the "linear" branch. 

Illustrations of how phase diagrams are modified as compared to the symmetric potential case are given in Fig. [SI 
E and F. As we increase rj from zero, a new double fixed point suddenly appears and bifurcates as we increase rj (see 
the next section). Thus we have two new fixed points (nolonger a pitchfork bifurcation). 



V. A JOSEPHSON JUNCTION APPROACH 



Now allow / to be time dependent and satisfy the one dimensional equation 



. df(x,t) 
1 dt 



-^ + V(x)+ V \f(x,t)f 



f(x,t), 



(27) 
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with potential V(x) in the form of a double well, which is not necessarily symmetric. To establish a link between 
the above equation and the Josephson model, we first focus on the energy spectrum of the system in the linear limit 
(r/ = 0). It consists of pairs of energy levels separated by a gap that is proportional to the height of the barrier; 
for a sufficiently high barrier the spacing between the pairs is larger than the spacing within the first pair. In this 
case we can construct a variational analysis based on the lowest pair of levels, and ip2- We assume that f(x, t) is 
normalized to unity and approximate it by f{x,t) ~ aL(t)wL(x) + aR(t)wn(x), where wl,r{x) are defined as 



1 



The eigenstates tpi(x) and ip2(%) are orthonormal. Amplitudes must satisfy |<il| 2 
derivatives are approximately given by 



ia L>R (t) = E Q a L , R (t) - Ka R , L (t) + U L , R \a L . R (t)\ 2 a L , R (t), 



(28) 

\a R \ 2 = 1 and their time 
(29) 



where 



,r( x ) 



dx 2 



V(x) 



wl,r{x), K 



w* L (x) 



dx 2 



+ V(x) 



w R (x)dx, U LfR = r) / \w LtR (x)\ 



Note that Eq is the common value of two expressions (for wl(x) and w R (x)). The Josephson equations, 



Az- 



(30) 
(31) 



will follow upon defining 



a-L,R = 



IT* 



eKp(i6 LtR ), cj> = L - Or, 



A = (U R + U L )/4K, A = (Ur — U L )/4K, 

and rescaling the time 2Kt — > t. Here A is the ratio of nonlinear coupling to tunneling and A is the difference in the 
depth of the wells. With our simplifications and substitutions the suitably normalized Hamiltonian of the system can 
be obtained in the form (see also 0)01) 



H/K = Eq/K — v 7 ! - z 2 cos. 



A , o 
2 (1 + Z 



Az, 



(32) 



The parameter A is positive for repulsive interaction {rj > 0) and negative for attractive interaction (j] < 0). Note 
the two symmetries: A — ► —A, cf> — ► <fi + tt, z — > —z and A — > —A, z — > — z, <f) — > —<f>. The first of these 
symmetries implies that completely solving for rj > gives the solution for rj < 0. 

These equations differ from those governing Josephsonian oscillations in superconducting junctions by two additional 
terms: one proportional to A which derives from the nonlinear interaction (it has the same sign as rf) and the constant 
A, owing its existence to the asymmetry of the potential. 

Consider the stationary solutions of the Josephson equations. From Eq. i|30|) we see that (j> — 0; ±7T. In Fig. |3| we 
see how to find them graphically. For A = there are always two solutions with z = 0. The other two solutions 
appear for nonzero z when A > 1 (or A < — 1). In the case of nonzero A there are also always at least two solutions. 
The other two appear above A equal A c = (1 — A 2 / 3 ) 3 / 2 . 

Having the stationary solutions we can draw the energy dependence on A, shown in Fig. The two lowest 
eigenvectors of the nonlinear Schrodinger equation (solid lines) are compared with those resulting from the Josephson 
Junction approach (dashed lines) for the case of equally deep wells. Notice the good agreement. 

Now consider the dynamics of the A = case. Constant energy contours in z, <p phase space followed by the system 
for positive A are shown in Fig. [SJ (A-D), one each for A < 1, 1 < A < 2 and two for A > 2, each of which is generic. 
The difference between the second and third case concerns the possibility of self trapping solutions oscillating about 
an average z such that 4> covers all possible values in the third case. However, the fixed point dynamics is common to 
the latter three cases. Fixed points are at: (1) z = 0, <j> = 0; (2) z = 0, <j> = ±7r; (3) z = yl — A -2 , tfi — i 71 ", A > 1. 
The latter pair bifurcate from the second point as we increase A through A = 1, see Fig. 3 for an illustration of how 
this happens. 

We will now look at the stability of the three classes of fixed points. Assume perturbations such that z^z + 5ze xt 
and 4> + S(pe xt . Simple calculations give values of A for the three categories: 
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FIG. 3: (Color online) Stationary solutions of the Josephson equations, represented by the intersections of the solid lines 
(y — zfc , z ) with the dashed ones (y = Az for A = 0.5; 1.8 - symmetric potential case) and dashed dotted line (y — Az + A 

for A = 0.5; 1.8 and A = 0.3- 



E 




FIG. 4: (Color online) Bifurcation diagram. The two sets of curves that almost coincide are obtained from the Josephson 
Hamiltonian (dashed lines) and the exact formulas of sections II and III (solid lines). Beyond the bifurcation point the 
symmetric solutions are unstable and the asymmetric ones are stable. Each point on the stable bifurcated branch corresponds 
to two mirror image solutions. 

(1) A 2 = —(1 + A) (phase point in the (<f>,z) plane moves on an ellipse like trajectory around (0,0)) 

(2) A 2 = A — 1 (fixed point stable when A < 1, but when A > 1, any perturbation moves out along an arm of the 
separatrix emerging from (±7r,0)) 

(3) A 2 = 1 — A 2 (phase point moves on an ellipse like trajectory around one of (±7t,±V1 — A -2 )) 

Thus, according to this criterion the first fixed point is always stable for A > — 1. The antisymmetric solution (2) is 
stable for A < 1 and unstable for A > 1. The bifurcated pair (3) is always stable and we have a typical pitchfork 
bifurcation at A = 1. These results are in full agreement with a numerical stability analysis based on the nonlinear 
Schrodinger equation (see Fig. El. We might add that they contradict some statements in the literature, e.g. 01 an d 

m 

One might wonder how the Josephson bifurcation picture ties up with the exact solutions considered earlier. Suppose 
we have an analytic solution given by mi, mz, ma, fi and d. As <f> = ±7r in our considerations, this solution clearly 
corresponds to one of the bifurcated fixed points in the Josephson model. How can we determine the corresponding 
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FIG. 5: Phase space diagrams. The first four frames correspond to symmetric potential wells, the latter two to a deeper well 
on the left (A = 0.3). Note the differences in the trajectories between the cases B and E, especially in the "waves" that cover 
the entire 0-range. The fixed points are still present in the corners of the fourth frame (D) but do not turn up on this scale. 



a) b) 




0.5 1 1.5 2 0.5 1 1.5 2 

t t 



c) d) 




FIG. 6: Evolution of z(t) for the antisymmetric states, (j> — n: a) r\ — —10 , A = —4.52 ; b) f] = — 1, A = —0.452 ; c) r\ = 1, 
A = 0.452; d) r\ = 10, A = 4.52. As we see in the case of a), b) and c) the solution is stable. The periods of oscillations match 
the formulas derived from the values of A given in Section 3, T = 27r/|A|. In case d) the solution is unstable. The stability of 
the asymmetric bifurcated branch has also been confirmed. 



value of A so that z takes the proper value? A good approximation to z when d is small, is (Appendix): 

_ [Qi - kl)oj - fci£(fciu;|mi) + fc 3 £'(fc 3 u;|m3)] 
Z ~ \[k\ + - kiE(kiu\mi) - k 3 E(k 3 uj\m 3 )] ^ ' 

and A = (1 — z 2 ) -1 / 2 . The apparent independence of d is deceptive, as d will determine mi and m 3 , above. In 
particular when d = 0, mi = m 3 , k% = fc 3 and z — as expected. The antisymmetric solution is recovered. 
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When A ^ 0, the fixed points are still located at (f> = 0; ±tt. The stationary values of z are now roots of the quartic 



AV + 2AAz 3 + (A 2 + 1 - A 2 )z 2 - 2AAz - A 2 = (34) 

and are shifted down as compared to the case of A = 0. Illustrations of how phase diagrams are modified as compared 
to the symmetric potential case are given in Fig. [3 E and F. Now phase curves covering all possible <fi values and 
such that z changes sign are possible. As we increase A from zero, a new double fixed point suddenly appears at a 
critical value of A equal A c = (1 — A 2 / 3 ) 3 / 2 , and bifurcates as we increase A. Thus for A > A c we have two new fixed 
points. A stability analysis yields results similar to the above, for A = 0, but values of \ are now given in terms of 
roots of the quartic Zi. 



VI. CONCLUSIONS 



In this paper we have thoroughly investigated the behavior of Bose-Einstein condensates in double square well 
potentials, both symmetric and asymmetric. A simple method for obtaining exact solutions for repulsive interaction 
was outlined (similarly as in [2] for attractive interaction). We treat the system both exactly and by a Josephson 
Junction model. We have checked the Josephson model results, both static and dynamic, against exact calculations. 
Agreement is suprisingly good. Some controversies about the stability, to be found in the literature, have been 
resolved. 



VII. ACKNOWLEDGEMENT 



The authors would like to acknowledge support from KBN Grant 2P03 B4325 (E.I.), KBN Grant 1P03B14629 (P.Z.) 
and the Polish Ministry of Scientific Research and Information Technology under grant PBZ-MIN-008/P03/2003 
(M.T.). 

Consultations with Professor George Rowlands were very helpful. 



VIII. APPENDIX 



A. Symmetry preserving case (symmetric wells) 



To find the A/i correction we eliminate A and A 2 from equations (JSJ and 10 

p _ kcn(kuj\m)dn(kLu\m) ^ k 2 dn(k 2 b\m 2 ) _ ^ 
sn(koj\m) sn(fc2&|m2)cn(fc26|m2) 

and calculate the perfect differential of F(m, m2,/i) for small 77. As the first two differentials follow from Eqs. I|14|) 
A/j, can be so obtained. This calculation is somewhat less straightforward. It completes the calculation of to, m 2 and 
fj, in the small 77 limit, our starting point. 

In the limit to and 1 — TO2 tending to zero Eq. (|15|l is recovered from l|35|) . The general equation for small increments 
of m, m2 and /1 is 

. OF OF OF 

AF = —Am + - — Am 2 + t^A/j = 
am am 2 u/x 

Am — to, Ato2 = to 2 — 1, Aji = i-i — /icb (36) 

and so 

/ OF 8F „ A fdFY 1 

^ = {-^ m+ ^ 1 - m ^)[^) ' (37) 

where in the perturbation limit to and 1 — TO2 are proportional to r\ and are given by equation (|14|l . We find after 



some calculations using known identities [T^j 



9F , /3„ 1 . „, 

— — = — k — 6 + — sin 2kuj 
dm \4 4 

^ 1 -5-J-L = - A " 



9/i 2fc~ 2fc 2 2 sin 2 few 

and 

6 = cot (few) — 



sin (Jeio) 

L = coth(fc 2 6) j . 

sinh (k 2 b) 



where k and k 2 are taken in the linear limit. 



B. Symmetry breaking case (symmetric wells) 

We obtain 

2 _ 2mife 2 2 _ 2m 3 fcg 2 _ 2(1 - m 2 )fcf 

A l — 7 ^3 — > A 2 — 1 

H = (1 + mi)fc 2 = (1 + m 3 )fc 2 = (m 2 - 2)fc 2 + V . 
The continuity conditions at x — ±6 are now generalized to: 

.9o = kiy/mi sn(fcxw|mi) — fc 2 -\/ (1 — m 2 ) sc(fc 2 (6 — rf)|m 2 ) = 

.9i = fc3\A«3 sn(fc 3 o;|m3) - fc 2V / (1 - m 2 ) sc(fc 2 (6 + d)|m 2 ) = 

.92 = fc 2 \/ TO i cn(feiw|mi) dn(fciw|r7ii) + fc 2 \/I — m 2 dc(fc 2 (fc — d)|m 2 ) nc(fc 2 (fo — rf)|m 2 ) 

.93 = fc 2 -\A™3 cn(/c 3 o;|?7i3) dn(fc 3 o;|r7i3) + fc 2 \/l — m 2 dc(fc 2 (6 + d)|m 2 ) nc(fc 2 (6 + d)|m 2 ) 

and the normalization condition is now: 

.94 = 2fc 2 [sn(fc 2 (6 - d)\m 2 ) dc(/c 2 (6 - d)\m 2 ) + sn(fc 2 (6 + d)\m 2 ) dc(fc 2 (6 + d)\m 2 )] 
-2k 2 [E{k 2 {b + d)\m 2 ) + E{k 2 (b - d)\m 2 )] 
+2k\ [kioj — E{k\ui\mi)] + 2k 3 [k 3 ui — E(k3U>\ms)] = 

If we can assume k 2 d much smaller than one, Eqs. I|41|) and l|42|l up to second order simplify to: 

(2 - m 2 )sc(fc 2 6|m 2 ) + 2(1 - m 2 )sc 3 (/c 2 ^|m 2 ) 

X(mi) - x(m 3 ) = 2D 2 rf 

m 2 

2D 

Xim-i) + x( m 3) = ; — dc(/c 2 6|m 2 )nc(fc 2 &|m 2 ) 

k 2 

ip(mi) — ip{m 3 ) — —2dk\\J\ — m 2 dc(fc 2 &|m2)nc(fc 2 6|m 2 ) 
4>(mi) + ip{m 3 ) — 2fc 2 \A — m 2 sc(k 2 b\m 2 ) 



T] = [<f>(mi) + 4>(m 3 )] + 4k 2 [sn(/c 2 6|m 2 )dc(/c 2 6|TO 2 ) - E(k 2 b\m 2 )] 
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where: 



X(m) 
tp(m) 



1 + m 



mcn(kui\m) dn(fcw|m) 



/i 



1 + m 
V -v 



m sn(fcw|ra) 



3/2 



2 - m 2 
2k [ku - E(ku\m)] 



\f\ — m 2 k 2 — 



1 + m 



TO 2 



and z — [4>(mi) — </>(m 3 )] / [(f>(mi) + (j)(m 3 )]. By comparing d determined by the first and third equations we can 
reduce the system to just four equations for four unknowns, mi, m 2: m ?> and //,. 
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